Introduction

This lab provides a fairly brief introduction to using R for Bayesian analysis. The first section will cover converting some of the models we have looked at in previous labs to Bayesian models using rstan (an interface to the stan language) and INLA. We’ll then go on to look at using INLA for spatial regression models.

R has a large number of add-on packages for Bayesian analysis listed here. Some notable ones are:

  • coda: provides routines for post-processing results
  • mcmcplots: provides prettier visualizations of MCMC results
  • MCMCpack: provides mainly standard statistical models in a simpler format

Installing INLA

The INLA package is not available through the usual CRAN repositories. Instead, you’ll need to install this manually. The following code will install the current stable version of INLA:

install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)

Once this is installed, we can load the libraries that we’ll need for this lab

library(dplyr)
library(ggplot2)
library(lme4)
library(sf)
library(INLA)
library(spdep)
library(tmap)

Non-spatial models

Linear regression

gap <- read.csv("../datafiles/gapminderData5.csv")
gap$year2 <- gap$year - 1952
gap$gdpPercap2 <- log(gap$gdpPercap)
fit_lm <- lm(lifeExp ~ gdpPercap2, gap)
summary(fit_lm)
## 
## Call:
## lm(formula = lifeExp ~ gdpPercap2, data = gap)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.778  -4.204   1.212   4.658  19.285 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -9.1009     1.2277  -7.413 1.93e-13 ***
## gdpPercap2    8.4051     0.1488  56.500  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.62 on 1702 degrees of freedom
## Multiple R-squared:  0.6522, Adjusted R-squared:  0.652 
## F-statistic:  3192 on 1 and 1702 DF,  p-value: < 2.2e-16

Basic INLA model

inla_lm <- inla(lifeExp ~ gdpPercap2, 
                data = gap)
inla_lm <- inla(lifeExp ~ gdpPercap2, 
                data = gap,
                control.compute = list(dic = TRUE, waic = TRUE),
                control.predictor = list(compute = TRUE))
summary(inla_lm)
## 
## Call:
##    c("inla(formula = lifeExp ~ gdpPercap2, data = gap, control.compute = 
##    list(dic = TRUE, ", " waic = TRUE), control.predictor = list(compute = 
##    TRUE))" ) 
## Time used:
##     Pre = 4.79, Running = 0.549, Post = 0.53, Total = 5.87 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -9.099 1.228    -11.510   -9.099     -6.691 -9.099   0
## gdpPercap2   8.405 0.149      8.113    8.405      8.697  8.405   0
## 
## Model hyperparameters:
##                                          mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.017 0.001      0.016    0.017
##                                         0.975quant  mode
## Precision for the Gaussian observations      0.018 0.017
## 
## Expected number of effective parameters(stdev): 2.00(0.00)
## Number of equivalent replicates : 851.93 
## 
## Deviance Information Criterion (DIC) ...............: 11760.48
## Deviance Information Criterion (DIC, saturated) ....: 1707.78
## Effective number of parameters .....................: 3.03
## 
## Watanabe-Akaike information criterion (WAIC) ...: 11760.97
## Effective number of parameters .................: 3.51
## 
## Marginal log-Likelihood:  -5899.79 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

In the linear model, we got an estimate of the residual standard deviation of around 7.6. This is a measure of the size of the errors from the model. In the INLA output, we get an estimate of the residual precision instead, which is the inverse of the variance. To compare, we can invert this and take the squre root to convert back to standard deviation units, and to check that the models are estimating the same value:

sqrt(1 / 0.017)
## [1] 7.66965
plot(inla_lm, plot.fixed.effects = TRUE,
     plot.random.effects = FALSE,
     plot.hyperparameters = FALSE,
     plot.predictor = FALSE)

plot(inla_lm, plot.fixed.effects = FALSE,
     plot.random.effects = FALSE,
     plot.hyperparameters = TRUE,
     plot.predictor = FALSE)

plot(gap$lifeExp, inla_lm$summary.linear.predictor$mean)
abline(0,1)

Generalized models

## Binomial regression
irished <- read.csv("../datafiles/irished.csv")
irished$DVRT.cen <- irished$DVRT - mean(irished$DVRT)
fit_glm <- glm(lvcert ~ DVRT.cen, 
               data = irished, 
               family = binomial(link = 'logit'))

summary(fit_glm)
## 
## Call:
## glm(formula = lvcert ~ DVRT.cen, family = binomial(link = "logit"), 
##     data = irished)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9077  -0.9810  -0.4543   1.0307   2.1552  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.278422   0.099665  -2.794  0.00521 ** 
## DVRT.cen     0.064369   0.007576   8.496  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 686.86  on 499  degrees of freedom
## Residual deviance: 593.77  on 498  degrees of freedom
## AIC: 597.77
## 
## Number of Fisher Scoring iterations: 3

Basic INLA model

inla_glm <- inla(lvcert ~ DVRT.cen, 
                 data = irished, 
                 family = "binomial",
                 control.compute = list(dic = TRUE, waic = TRUE),
                 control.predictor = list(compute = TRUE))
summary(inla_glm)
## 
## Call:
##    c("inla(formula = lvcert ~ DVRT.cen, family = \"binomial\", data = 
##    irished, ", " control.compute = list(dic = TRUE, waic = TRUE), 
##    control.predictor = list(compute = TRUE))" ) 
## Time used:
##     Pre = 5.54, Running = 0.323, Post = 0.383, Total = 6.24 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.279 0.100     -0.476   -0.279     -0.084 -0.278   0
## DVRT.cen     0.065 0.008      0.050    0.064      0.080  0.064   0
## 
## Expected number of effective parameters(stdev): 2.00(0.00)
## Number of equivalent replicates : 249.92 
## 
## Deviance Information Criterion (DIC) ...............: 597.75
## Deviance Information Criterion (DIC, saturated) ....: 597.75
## Effective number of parameters .....................: 1.99
## 
## Watanabe-Akaike information criterion (WAIC) ...: 597.79
## Effective number of parameters .................: 2.02
## 
## Marginal log-Likelihood:  -306.61 
## Posterior marginals for the linear predictor and
##  the fitted values are computed
plot(inla_glm, plot.fixed.effects = TRUE,
     plot.random.effects = FALSE,
     plot.hyperparameters = FALSE,
     plot.predictor = FALSE)

Mixed-effects models

fit_lmer <- lmer(lifeExp ~ year2 + (1 | country), gap)
summary(fit_lmer)
## Linear mixed model fit by REML ['lmerMod']
## Formula: lifeExp ~ year2 + (1 | country)
##    Data: gap
## 
## REML criterion at convergence: 9866.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.1692 -0.4910  0.1216  0.6030  2.4194 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  country  (Intercept) 123.15   11.097  
##  Residual              12.85    3.584  
## Number of obs: 1704, groups:  country, 142
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 50.51208    0.94547   53.43
## year2        0.32590    0.00503   64.79
## 
## Correlation of Fixed Effects:
##       (Intr)
## year2 -0.146
inla_lmer <- inla(lifeExp ~ year2 + f(country, model = "iid"), 
                  data = gap)
summary(inla_lmer)
## 
## Call:
##    "inla(formula = lifeExp ~ year2 + f(country, model = \"iid\"), data = 
##    gap)" 
## Time used:
##     Pre = 6.04, Running = 0.558, Post = 0.343, Total = 6.94 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) 50.512 0.947     48.651   50.512     52.372 50.512   0
## year2        0.326 0.005      0.316    0.326      0.336  0.326   0
## 
## Random effects:
##   Name     Model
##     country IID model
## 
## Model hyperparameters:
##                                          mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.078 0.003      0.072    0.078
## Precision for country                   0.008 0.001      0.006    0.008
##                                         0.975quant  mode
## Precision for the Gaussian observations      0.083 0.078
## Precision for country                        0.010 0.008
## 
## Expected number of effective parameters(stdev): 141.77(0.154)
## Number of equivalent replicates : 12.02 
## 
## Marginal log-Likelihood:  -4968.37
plot(inla_lmer, plot.fixed.effects = TRUE,
     plot.random.effects = FALSE,
     plot.hyperparameters = FALSE,
     plot.predictor = FALSE)

plot(inla_lmer, plot.fixed.effects = FALSE,
     plot.random.effects = TRUE,
     plot.hyperparameters = FALSE,
     plot.predictor = FALSE, cex = 1)

plot(inla_lmer, plot.fixed.effects = FALSE,
     plot.random.effects = FALSE,
     plot.hyperparameters = TRUE,
     plot.predictor = FALSE, 
     plot.prior = TRUE)

Spatial models

Read in the data

ufos <- st_read("../datafiles/nuforc/ufo_sf_annual.shp")
## Reading layer `ufo_sf_annual' from data source 
##   `/Users/u0784726/Dropbox/Data/devtools/geog6000/datafiles/nuforc/ufo_sf_annual.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1421 features and 44 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.71 ymin: 24.54235 xmax: -66.98702 ymax: 49.36967
## Geodetic CRS:  WGS 84

Select only 2014

ufos_2014 <- ufos %>%
  filter(year == 2014)
p1 <- tm_shape(ufos_2014) + tm_fill("count", style = "jenks") +
  tm_borders() +
  tm_layout(main.title = "UFO sightings 2014")

p2 <- tm_shape(ufos_2014) + tm_fill("population", style = "jenks") +
  tm_borders() +
  tm_layout(main.title = "State population 2014")

tmap_arrange(p1, p2)

Basic linear regression

ufo_lm <- lm(count ~ population, ufos_2014)
summary(ufo_lm)
## 
## Call:
## lm(formula = count ~ population, data = ufos_2014)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -261.79  -30.40  -12.68   22.04  219.68 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.524562  14.802688   2.332    0.024 *  
## population   0.019879   0.001543  12.882   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 76.66 on 47 degrees of freedom
## Multiple R-squared:  0.7793, Adjusted R-squared:  0.7746 
## F-statistic: 165.9 on 1 and 47 DF,  p-value: < 2.2e-16
ufo_inla_lm <- inla(count ~ population,
                    data = ufos_2014,
                    control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE),
                    control.predictor = list(compute = TRUE))
summary(ufo_inla_lm)
## 
## Call:
##    c("inla(formula = count ~ population, data = ufos_2014, control.compute 
##    = list(dic = TRUE, ", " waic = TRUE, cpo = TRUE), control.predictor = 
##    list(compute = TRUE))" ) 
## Time used:
##     Pre = 5.03, Running = 0.313, Post = 0.351, Total = 5.7 
## Fixed effects:
##               mean     sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) 34.525 14.936      5.058   34.524     63.964 34.525   0
## population   0.020  0.002      0.017    0.020      0.023  0.020   0
## 
## Model hyperparameters:
##                                         mean   sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.00 0.00       0.00     0.00
##                                         0.975quant mode
## Precision for the Gaussian observations       0.00 0.00
## 
## Expected number of effective parameters(stdev): 2.00(0.00)
## Number of equivalent replicates : 24.50 
## 
## Deviance Information Criterion (DIC) ...............: 568.69
## Deviance Information Criterion (DIC, saturated) ....: 58.39
## Effective number of parameters .....................: 3.16
## 
## Watanabe-Akaike information criterion (WAIC) ...: 574.64
## Effective number of parameters .................: 7.59
## 
## Marginal log-Likelihood:  -307.59 
## CPO and PIT are computed
## 
## Posterior marginals for the linear predictor and
##  the fitted values are computed
ufo_inla_lm$summary.fixed
##                    mean           sd 0.025quant    0.5quant  0.975quant
## (Intercept) 34.52460098 14.935909271 5.05803232 34.52414369 63.96379416
## population   0.01987854  0.001557171 0.01680645  0.01987849  0.02294778
##                    mode          kld
## (Intercept) 34.52460098 1.497393e-05
## population   0.01987854 1.500443e-05

Posterior distribution

ufo_inla_lm$marginals.fixed
ufo_inla_lm$marginals.hyperpar

Plot posterior distribution of intercept

plot_df <- as.data.frame(inla.tmarginal(function(x) x, 
                                        ufo_inla_lm$marginals.fixed$`(Intercept)`))
ggplot(plot_df, aes(x = x, y = y)) +
  geom_line(size = 2) +
  scale_x_continuous("Intercept") +
  scale_y_continuous("P") +
  theme_bw()

Plot posterior distribution of population coefficient

plot_df <- as.data.frame(inla.tmarginal(function(x) x, 
                                        ufo_inla_lm$marginals.fixed$population))
ggplot(plot_df, aes(x = x, y = y)) +
  geom_line(size = 2) +
  scale_x_continuous("Population") +
  scale_y_continuous("P") +
  theme_bw()

Plot posterior distribution of residual variance

plot_df <- as.data.frame(inla.tmarginal(function(x) 1/x, 
                                        ufo_inla_lm$marginals.hyperpar$`Precision for the Gaussian observations`))
ggplot(plot_df, aes(x = x, y = y)) +
  geom_line(size = 2) +
  scale_x_continuous("Residual variance") +
  scale_y_continuous("P") +
  theme_bw()

Add an informed prior on slope

ufo_inla_lm2 <- inla(count ~ population,
                     data = ufos_2014,
                     control.fixed = list(mean = 0, prec = 0.2, 
                                          mean.intercept = 0, prec.intercept = 0.2),
                     control.family=list(hyper=list(prec=list(prior='loggamma', 
                                                              param=c(0.1,0.1)))),
                     control.compute = list(dic = TRUE, waic = TRUE),
                     control.predictor = list(compute = TRUE))
summary(ufo_inla_lm2)
## 
## Call:
##    c("inla(formula = count ~ population, data = ufos_2014, control.compute 
##    = list(dic = TRUE, ", " waic = TRUE), control.predictor = list(compute 
##    = TRUE), control.family = list(hyper = list(prec = list(prior = 
##    \"loggamma\", ", " param = c(0.1, 0.1)))), control.fixed = list(mean = 
##    0, prec = 0.2, ", " mean.intercept = 0, prec.intercept = 0.2))") 
## Time used:
##     Pre = 5.03, Running = 0.259, Post = 0.346, Total = 5.63 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.713 2.217     -3.642    0.713      5.062 0.713   0
## population  0.022 0.001      0.020    0.022      0.025 0.022   0
## 
## Model hyperparameters:
##                                         mean   sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.00 0.00       0.00     0.00
##                                         0.975quant mode
## Precision for the Gaussian observations       0.00 0.00
## 
## Expected number of effective parameters(stdev): 1.02(0.004)
## Number of equivalent replicates : 48.01 
## 
## Deviance Information Criterion (DIC) ...............: 571.62
## Deviance Information Criterion (DIC, saturated) ....: 51.66
## Effective number of parameters .....................: 2.08
## 
## Watanabe-Akaike information criterion (WAIC) ...: 577.02
## Effective number of parameters .................: 6.33
## 
## Marginal log-Likelihood:  -295.48 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

Spatial neighborhood and weight matrix

ufos_2014$idarea <- 1:nrow(ufos_2014)
ufos_nb <- poly2nb(ufos_2014)
ufos_W <- nb2mat(ufos_nb)
plot(st_geometry(ufos_2014), reset = FALSE)
plot(ufos_nb, st_coordinates(st_centroid(ufos_2014)), add = TRUE, col = 2, lwd = 2)

ufo_inla_bym <- inla(count ~ population + f(idarea, model = "bym2", graph = ufos_W), 
                     data = ufos_2014,
                     control.compute = list(dic = TRUE, waic = TRUE),
                     control.predictor = list(compute = TRUE)
)
summary(ufo_inla_bym)
## 
## Call:
##    c("inla(formula = count ~ population + f(idarea, model = \"bym2\", ", " 
##    graph = ufos_W), data = ufos_2014, control.compute = list(dic = TRUE, 
##    ", " waic = TRUE), control.predictor = list(compute = TRUE))" ) 
## Time used:
##     Pre = 28.2, Running = 0.478, Post = 0.402, Total = 29.1 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) 40.013 6.129     27.891   40.041      51.97 40.097   0
## population   0.019 0.001      0.018    0.019       0.02  0.019   0
## 
## Random effects:
##   Name     Model
##     idarea BYM2 model
## 
## Model hyperparameters:
##                                             mean       sd 0.025quant 0.5quant
## Precision for the Gaussian observations 3.08e+04 1.48e+04   9292.433 2.86e+04
## Precision for idarea                    1.00e-03 0.00e+00      0.001 1.00e-03
## Phi for idarea                          3.54e-01 1.03e-01      0.188 3.41e-01
##                                         0.975quant     mode
## Precision for the Gaussian observations   6.56e+04 2.28e+04
## Precision for idarea                      1.00e-03 1.00e-03
## Phi for idarea                            5.84e-01 3.05e-01
## 
## Expected number of effective parameters(stdev): 49.00(0.00)
## Number of equivalent replicates : 1.00 
## 
## Deviance Information Criterion (DIC) ...............: -295.17
## Deviance Information Criterion (DIC, saturated) ....: 125.63
## Effective number of parameters .....................: 62.82
## 
## Watanabe-Akaike information criterion (WAIC) ...: -318.00
## Effective number of parameters .................: 30.41
## 
## Marginal log-Likelihood:  -463.53 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

Add a prior on the spatial term.

prior_bym2 <- list(
  prec = list(
    prior = "pc.prec",
    param = c(1, 0.01)),
  phi = list(
    prior = "pc",
    param = c(0.5, 2 / 3))
)
ufo_inla_bym2 <- inla(count ~ population + f(idarea, model = "bym2", 
                                            graph = ufos_W,
                                            hyper = prior_bym2), 
                     data = ufos_2014,
                     control.compute = list(dic = TRUE, waic = TRUE),
                     control.predictor = list(compute = TRUE),
                     control.inla = list(strategy = "adaptive",int.strategy = "eb")
)
summary(ufo_inla_bym2)
## 
## Call:
##    c("inla(formula = count ~ population + f(idarea, model = \"bym2\", ", " 
##    graph = ufos_W, hyper = prior_bym2), data = ufos_2014, control.compute 
##    = list(dic = TRUE, ", " waic = TRUE), control.predictor = list(compute 
##    = TRUE), control.inla = list(strategy = \"adaptive\", ", " int.strategy 
##    = \"eb\"))") 
## Time used:
##     Pre = 25.8, Running = 0.461, Post = 0.368, Total = 26.6 
## Fixed effects:
##               mean     sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) 34.525 14.124      6.794   34.524     62.233 34.525   0
## population   0.020  0.001      0.017    0.020      0.023  0.020   0
## 
## Random effects:
##   Name     Model
##     idarea BYM2 model
## 
## Model hyperparameters:
##                                          mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.000 0.000      0.000    0.000
## Precision for idarea                    1.192 0.413      0.587    1.124
## Phi for idarea                          0.048 0.016      0.022    0.047
##                                         0.975quant  mode
## Precision for the Gaussian observations      0.000 0.000
## Precision for idarea                         2.187 1.001
## Phi for idarea                               0.082 0.044
## 
## Expected number of effective parameters(stdev): 2.01(0.00)
## Number of equivalent replicates : 24.40 
## 
## Deviance Information Criterion (DIC) ...............: 566.34
## Deviance Information Criterion (DIC, saturated) ....: 55.66
## Effective number of parameters .....................: 2.01
## 
## Watanabe-Akaike information criterion (WAIC) ...: 569.79
## Effective number of parameters .................: 4.78
## 
## Marginal log-Likelihood:  -283.33 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

Plot spatial effect

ufos_2014$u <- ufo_inla_bym2$summary.random$idarea$mean[1:49]
tm_shape(ufos_2014) + tm_fill("u", style = "jenks") +
  tm_borders() +
  tm_layout(main.title = "State population 2014")

Plot Bayesian estimates of sightings

ufos_2014$yhat <- ufo_inla_bym2$summary.linear.predictor$`0.025quant`
tm_shape(ufos_2014) + tm_fill("yhat", style = "jenks") +
  tm_borders() +
  tm_layout(main.title = "State population 2014")

Rate model

ufos_2014$Ei <- (sum(ufos_2014$count) / sum(ufos_2014$population)) * ufos_2014$population
ufos_2014$SIR <- ufos_2014$count / ufos_2014$Ei
tm_shape(ufos_2014) + tm_fill("SIR") + 
  tm_borders() + tm_layout(main.title = "UFO SIR 2014")

ufo_inla_rr <- inla(count ~ 1 + f(idarea, model = "bym2", 
                           graph = ufos_W,
                           hyper = prior_bym2), 
             data = ufos_2014, family = "poisson", E = Ei, 
             control.compute = list(dic = TRUE, waic = TRUE),
             control.predictor = list(compute = TRUE),
             control.inla = list(strategy = "adaptive",int.strategy = "eb")
)
summary(ufo_inla_rr)
## 
## Call:
##    c("inla(formula = count ~ 1 + f(idarea, model = \"bym2\", graph = 
##    ufos_W, ", " hyper = prior_bym2), family = \"poisson\", data = 
##    ufos_2014, ", " E = Ei, control.compute = list(dic = TRUE, waic = 
##    TRUE), ", " control.predictor = list(compute = TRUE), control.inla = 
##    list(strategy = \"adaptive\", ", " int.strategy = \"eb\"))") 
## Time used:
##     Pre = 25.6, Running = 0.379, Post = 0.36, Total = 26.3 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.033 0.047      -0.06    0.033      0.126 0.033   0
## 
## Random effects:
##   Name     Model
##     idarea BYM2 model
## 
## Model hyperparameters:
##                       mean   sd 0.025quant 0.5quant 0.975quant  mode
## Precision for idarea 4.759 1.30      2.660    4.610      7.738 4.328
## Phi for idarea       0.568 0.24      0.115    0.584      0.949 0.771
## 
## Expected number of effective parameters(stdev): 45.04(0.00)
## Number of equivalent replicates : 1.09 
## 
## Deviance Information Criterion (DIC) ...............: 420.40
## Deviance Information Criterion (DIC, saturated) ....: 106.46
## Effective number of parameters .....................: 45.18
## 
## Watanabe-Akaike information criterion (WAIC) ...: 414.14
## Effective number of parameters .................: 28.14
## 
## Marginal log-Likelihood:  -236.32 
## Posterior marginals for the linear predictor and
##  the fitted values are computed
ufos_2014$RR <- ufo_inla_rr$summary.fitted.values[, "mean"]

tm_shape(ufos_2014) + tm_fill("RR", palette = "-inferno") +
  tm_borders() +
  tm_layout(main.title = "RR")

ufos_2014$exc <- sapply(ufo_inla_rr$marginals.fitted.values,
                        function(marg){1 - inla.pmarginal(q = 2, marginal = marg)})

tm_shape(ufos_2014) + tm_fill("exc", palette = "Blues") +
  tm_borders() +
  tm_layout(main.title = "Excedence probability (RR > 2)")

airports <- read.csv("../datafiles/nuforc/airports.csv")
airports$lairports <- log(airports$airports)
ufos_2014 <- merge( ufos_2014, airports, by.x = "postal", by.y = "state")
ggplot(ufos_2014, aes(x = airports, y = count)) +
  scale_x_log10("Airports") + scale_y_log10("Sightings") +
  geom_point()

prior_bym2 <- list(
  prec = list(
    prior = "pc.prec",
    param = c(0.6 / 0.31, 0.01)),
  phi = list(
    prior = "pc",
    param = c(0.005, 1 / 10))
)
ufo_inla_rr2 <- inla(count ~ lairports + 
               f(idarea, model = "bym2", 
                 graph = ufos_W,
                 hyper = prior_bym2), 
             data = ufos_2014, family = "poisson", E = Ei, 
             control.compute = list(dic = TRUE, waic = TRUE),
             control.predictor = list(compute = TRUE),
             control.inla = list(strategy = "adaptive",int.strategy = "eb")
)
summary(ufo_inla_rr2)
## 
## Call:
##    c("inla(formula = count ~ lairports + f(idarea, model = \"bym2\", ", " 
##    graph = ufos_W, hyper = prior_bym2), family = \"poisson\", ", " data = 
##    ufos_2014, E = Ei, control.compute = list(dic = TRUE, ", " waic = 
##    TRUE), control.predictor = list(compute = TRUE), ", " control.inla = 
##    list(strategy = \"adaptive\", int.strategy = \"eb\"))" ) 
## Time used:
##     Pre = 23.2, Running = 0.344, Post = 0.362, Total = 23.9 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)  0.419 0.482     -0.529    0.420      1.364  0.421   0
## lairports   -0.068 0.084     -0.233   -0.068      0.097 -0.068   0
## 
## Random effects:
##   Name     Model
##     idarea BYM2 model
## 
## Model hyperparameters:
##                       mean    sd 0.025quant 0.5quant 0.975quant mode
## Precision for idarea 4.915 1.212      2.915    4.792      7.638 4.56
## Phi for idarea       0.142 0.144      0.005    0.093      0.542 0.01
## 
## Expected number of effective parameters(stdev): 45.72(0.00)
## Number of equivalent replicates : 1.07 
## 
## Deviance Information Criterion (DIC) ...............: 420.62
## Deviance Information Criterion (DIC, saturated) ....: 106.67
## Effective number of parameters .....................: 45.87
## 
## Watanabe-Akaike information criterion (WAIC) ...: 414.30
## Effective number of parameters .................: 28.51
## 
## Marginal log-Likelihood:  -242.86 
## Posterior marginals for the linear predictor and
##  the fitted values are computed